
## This code estimates the Asahi Shimbun data using ordinal variational IRT
## Takes dat_full.csv and prep.analysis.R as input files
## Outputs "out.VI3cat.1dim.PV0.RData" file, which is used for graphics in "figures.R"
## Author: Yuki Shiraito, 9/29/2014

dat.all <- read.csv("dat_full.csv")[, -1]

## Originally functions in separate packages, now all calls should be to emIRT
#library(ordIRT)
library(emIRT)

### number of dimensions; 1L or 2L
n.dims <- 1L

### data to be used; 0: both, 1: politicians only, 2: voters only
pv <- 0L

### preprocess data
source("prep.analysis.R", echo=TRUE)

### drop items with 2 response categories
drop.items <- apply(dat.all, 2, function(x){length(sort(unique(x)))}) == 2
dat.all <- dat.all[, !drop.items]

### convert the responses to 3 categories
dat.all <- apply(dat.all, 2, function(x) {
    if (length(sort(unique(x))) == 4) {
        temp.logical1 <- x == sort(unique(x))[1]
        temp.logical2 <- x == sort(unique(x))[2]
        temp.logical3 <- x %in% sort(unique(x))[3:4]
        x <- ifelse(temp.logical1, 1, x)
        x <- ifelse(temp.logical2, 2, x)
        x <- ifelse(temp.logical3, 3, x)
    } else if (length(sort(unique(x))) == 5) {
        temp.logical1 <- x %in% sort(unique(x))[1:2]
        temp.logical2 <- x == sort(unique(x))[3]
        temp.logical3 <- x %in% sort(unique(x))[4:5]
        x <- ifelse(temp.logical1, 1, x)
        x <- ifelse(temp.logical2, 2, x)
        x <- ifelse(temp.logical3, 3, x)
    }
})

### dimension of data
n.obs <- nrow(dat.all)
n.votes <- ncol(dat.all)

### generate priors and starting values
priors <- makePriors(n.obs, n.votes, n.dims)
start.values <- getStarts(n.obs, n.votes, n.dims)
start.values$DD <- matrix(rep(0.5,n.votes),ncol=1)
start.values$tau <- matrix(rep(-0.5,n.votes),ncol=1)
start.values$beta <- matrix(runif(n.votes, -1, 1), ncol = 1)


### Modify to eliminate alpha
priors$beta$mu <- matrix(0,1,1)
priors$beta$sigma <- matrix(25,1,1)

### "rollcall" object
dat.all <- as.matrix(dat.all)
dat.all[is.na(dat.all)] <- 0
## ### objects to be supplied to fastest()
## priors <- makePriors(n.obs, n.votes, n.dims)
## if (n.dims == 2) {
##     priors$x$sigma[2, 2] <- 10
## }
## start.values <- getStarts(n.obs, n.votes, n.dims)


## For new update
priors$beta$mu <- matrix(0,2,1)
priors$beta$sigma <- diag(25,2)

### variational inference
time <- system.time({
out.varinf <- ordIRT(.rc = dat.all, .starts = start.values, .priors = priors,
                     .D = n.dims, .control = {list(#threads = 8,
                         verbose = TRUE,
                         thresh = 1e-6, maxit = 2500)})
})


save.image(paste("out.VI3cat.", n.dims, "dim.PV", pv, ".RData", sep = ""))
